#include <iostream>
#include <vector>
#include <string>
#include <fstream>
#include <map>
#include <stdlib.h>
#include <lineages.h>
#include <algorithm>
#include <stdio.h>
#include <sstream>
#include <math.h>
#include <time.h>
#include <shared.h>

using namespace std;

int fitnessType;
int robustType;
int disturbType;

double baseDelRate = 1.0;

void usage()
{
	cerr << "Calculate the robustness, complexity, etc for a collection of lineages" << endl;
	cerr << "Usage: misc lin fit robust disturb start num seed collect" << endl;
	cerr << "Example: misc fun.alm 0 1 22 1 100 123 fun.random.lineages" << endl;
	cerr << "Parameters (all are needed) :"<< endl;
	cerr << " lin       : Name of a lineage file (one cell per line)" << endl;
	cerr << "             The first row of the output file will be the result for this lineage." << endl;
	cerr << "             The program use a silly match-by-string way to recognize what lineages are being used." << endl;
	cerr << "             The matches are as follows:" << endl;
	cerr << "             \"fun\".  C. elegans lineage and functional classification" << endl;
	cerr << "             \"pma\".  P. marina lineage and functional classification" << endl;
	cerr << "             \"hro\".  H. roretzi lineage and functional classification" << endl;
	cerr << "             \"cel\".  C. elegans lineage and functional classification, with neuron subtypes defined" << endl;
	cerr << "             \"cpe\".  C. elegans lineage up to the post-embryonic stage, and functional classification" << endl;
	cerr << "             So your filename must contain one of these five strings, otherwise the program will fail" << endl;
	cerr << "             The program need to know what lineage is being used because" << endl;
	cerr << "             it need to \"understand\" the names of functional classification" << endl;
	cerr << "             I recommend only those lineage files downloadable from the program's google code project site." << endl;
    cerr << " fit       : The type of \"fitness\" to be calculated." << endl;
    cerr << "             Could be 0-4" << endl;
	cerr << "             0.  robustness" << endl;
	cerr << "             1.  (1-complexity)" << endl;
	cerr << "             2.  robustness/complexity" << endl;
	cerr << "             3.  robustness^5/complexity" << endl;
	cerr << "             4.  robustness^100/complexity" << endl;
	cerr << " robust    : If you used 0 or 2-4 for \"fit\", this parameter determines the type of robustness to be used." << endl;
	cerr << "             Could be 1-2" << endl;
	cerr << "             1.  product of fractional survival in each type, ignoring apoptotic cells (used in our paper)" << endl;
	cerr << "             2.  fractional survival of all non-apoptotic cells" << endl;
	cerr << " disturb   : Type of disturbance happen on the lineage. The robustness against this disturbance will be used." << endl;
	cerr << "             Could be 1 or 22" << endl;
	cerr << "             1.  a single necrosis among non-root cells" << endl;
	cerr << "             22.  program failure at the rate of 1/N_i, N_i is the number of internal nodes" << endl;
	cerr << " start     : Start with this lineage within the collection" << endl;
	cerr << " num       : Deal with this many lineages in the collection" << endl;
	cerr << " seed      : Seed for randomization" << endl;
	cerr << " collect   : File name for the lineage collection (one lineage per line)" << endl;
	cerr << " " << endl;
	cerr << "Output would be on STDOUT" << endl;
	cerr << " " << endl;
	return;
}

int main(int argc, char *argv[])
{
	if(argc != 9)
	{
		usage();
		exit(255);
	}
    string realLineagesFile = argv[1];
    fitnessType = atoi(argv[2]);
    robustType = atoi(argv[3]);
    disturbType = atoi(argv[4]);
    int simStart = atoi(argv[5]);
    int simLoop = atoi(argv[6]);
    int randSeed = atoi(argv[7]);
    string simFile = argv[8];
    srand(randSeed);

    double realComplexity,realRobustness;
    double tmpDbl = 0;
    vector<int> termTypes;
    vector<string> termLineages;
    vector<vector<string> > lineagesPool;
    vector<vector<int> > typesPool;

    loadRealLineages(realLineagesFile,termTypes,termLineages);
    Lineages realLin(termTypes,termLineages);

    realLin.getDepthStat();
    realComplexity = realLin.complexity();
    realLin.ruleDeletionRate = baseDelRate / (realLin.termTypes.size() - 1);
    // Here the denominator is (realLin.termTypes.size() - 1)
    // This is correct!! See explanations below, when calculating this metric for the simulated trees.
    realRobustness = realLin.robustness(disturbType,robustType);

    if (disturbType == 2)
    {
        realRobustness = (1 - baseDelRate ) + baseDelRate * realRobustness;
	}
    cout << "Parameters: FitnessType(" << fitnessType;
    cout << ") ; RobustType(" << robustType << ")";
    cout << endl;

    cout << "Lineage Name\tComplexity\tRobustness\tMax depth\tMean depth\tS.D. depth\tfitness" << endl;
    cout << realLineagesFile << "\t" << realComplexity << "\t" << realRobustness << "\t" << realLin.depthMax << "\t";
    cout << realLin.depthMean << "\t" << realLin.depthSD << "\t" << fitness(fitnessType,realRobustness,realComplexity) << endl;

    for(int i=simStart; i<simStart+simLoop; i++)
    {
        lineagesPool.clear();
        typesPool.clear();
        loadLineageCollection(simFile,i,1,typesPool,lineagesPool);

        termTypes = typesPool[0];
        termLineages = lineagesPool[0];

        Lineages simLin(termTypes,termLineages);
        double simComplexity,simRobustness;
        simLin.getDepthStat();
        simComplexity = simLin.complexity();
        simLin.ruleDeletionRate = baseDelRate / (simLin.termTypes.size() - 1);
        // Here the denominator is (simLin.termTypes.size() - 1)
        // This is correct!! Don't use simLin.internalTypeCnt!!
        // Because rule deletion reflex genetic mutation rate. And in a lineage with more rules, it should have higher chances of mutation/rule deletion.
        // If simLin.internalTypeCnt is used here, a lineage with more rules with have lower mutation rate. That is not what we want.
        // On the other hand, this could of couse be set to some different value, as long as they are reasonable and the same across different simulated lineages based on the same real lineage (Stay the same in the group that with are comparing within) .
        simRobustness = simLin.robustness(disturbType,robustType);

        if(disturbType == 2) // This is obsolete. Don't use it.
        {
            tmpDbl = simLin.internalTypeCnt / realLin.internalTypeCnt ;
            if(baseDelRate * tmpDbl > 1)
            {
                cerr << "Basal Deletion Rate is too high for rule deletion" << endl;
                cerr << "Disturb type (2) is only intended for rule deletion rate < 1" << endl;
                cerr << "Here the effective rule deletion rate for this simulated lineages " << i << " is " << baseDelRate * tmpDbl << endl;
                cerr << "Please consider reduce it" << endl;
            }
            simRobustness = (1 - baseDelRate * tmpDbl) + baseDelRate * tmpDbl * simRobustness;
        }
        cout << i << "\t" << simComplexity << "\t" << simRobustness << "\t" << simLin.depthMax << "\t";
        cout << simLin.depthMean << "\t" << simLin.depthSD << "\t" << fitness(fitnessType,simRobustness,simComplexity)<< endl;
    }

    return 0;
}
